gusucode.com > 《MATLAB图像与视频处理实用案例详解》代码 > 《MATLAB图像与视频处理实用案例详解》代码/第 19 章 基于语音识别的信号灯图像模拟控制技术/voicebox/modspect.m
function [c,qq,ff,tt,po]=modspect(s,fs,m,nf,nq,p) %MODSPECT Calculate the modulation spectrum of a signal C=(S,FS,W,NC,P,N,INC,FL,FH) % % % Simple use: % % Inputs: % s speech signal % fs sample rate in Hz (default 11025) % m mode string containing any sensible combination of the following letters; these % options override any specified in argument p. % % g apply AGC to the input speech % % r/n/m time domain windows: rectangular/hanning/hamming [def: hamming] % T/N/M mel/mod-mel domain windows: triangular/hanning/hamming [def: triangular] % % a/c mel filters act in the amplitude or complex domain [def: power] % p/l mel outputs are in power or log power [def: amplitude] % A mod-mel filters act in the power domain [def: power] % P/L mod-mel outputs are in amplitude or log power [def: amplitude] % k include DC component of modulation spectrum [not implemented] % % F Take a DCT in the mel direction % Z include 0'th order mel-cepstral coef % % f Take a DCT in the mod-mel direction % z include 0'th order mod-mel cepstral coef % % D include d/df coefficients in * per mel (not valid with 'F') % d include d/dt coefficients in * per second % % u give frequency axes in mels and mod-mels rather than in Hz % % s plot the result % S plot intermediate results also % % nf four-element vector giving: % nf(1) = number of mel bins (or increment in k-mel) [0.1] % nf(2) = min filterbank frequency [40 Hz] % nf(3) = max filterbank frequency [10kHz] % nf(4) = number of DCT coefficientss excl 0'th [25]. % Omitted values use defaults. % nq four-element vector giving % nq(1) = number of mel-mod bins (or increment in 10 mod-mel units) % nq(2) = min filterbank frequency [50 Hz] % nq(3) = max filterbank frequency [400 Hz] % nq(4) = number of DCT coefficients excl 0'th [15] % Omitted values use defaults. Note that the mel-mod frequency scale % is like the mel scale but reduced by a factor of 100. % I.e. mod-mel(f)=mel(100 f)/100; thus 10Hz corresponds to 10 mod-mels. % p parameter structure containing some or all of the parameters listed % at the start of the program code below. p may be given as the % last input parameter even if some or all of the inputs m, nf and nq % have been omitted. % % % % Outputs: % c[*,nf,nt] modulation spectrum where nt the number of time frames and % nf is the normally number of mel frequency bins (but reduced by 1 if % 'F' is specified without 'Z'). % The middle index is for the output coefficients % which are concatenated in the order [base, delta-p, delta-t]. % The number of base coefficients is normally the % number of mon-mel filters but will be reduced by 1 % if 'f' is specified but 'z' is not. % qq centre frequencies of modulation bins before any DCT (Hz or mel if 'u' set) % ff centre frequencies of freqiency bins before any DCT (Hz or mel if 'u' set) % tt time axis (sample s(i) is at time i/fs) % po a copy of the parameter structure that was actually used. % % Algorithm Outline [parameters are in brackets]: % % (1) Apply AGC to the speech to preserve a constant RMS level [tagc,tagt;'g'] % (2) Divide into overlapping frames [tinc,tovl], window [twin;'rnm'], zero-pad [tpad,trnd] and FFT % (3) Apply mel filters in power or amplitude domain [fmin,fmax,fbin,fwin,fpow;'TNMac'=nf] % (4) Regarding each filterbank output as a time-varying signal power or amplitude [ppow;'pl'], % divide into overlapping frames [pinc,povl], window [pwin;'rnm'], zero-pad [ppad,prnd] and FFT % (5) Apply mod-mel filters in power or amplitude domain [mmin,mmax,mbin,mwin,mpow;'TNMA'=nq] % (6) Take optional log [qpow;'PL'] % (7) Optionally apply a DCT in the mel [qdcp;'F'] and/or mod-mel [qdcq;'f'] directions preserving % the only some low quefrency coefficients [ dncp, dzcp;'Z'=nf, dncq, dzcq;'z'=nq] % (8) Calculate derivatives over mel-frequency [ddep,ddnp;'D'] and/or over time [ddet,ddnt;'d'] % % Notes: All logs are limitied in dynamic range [logr], output units can be Hz or mel [unit;'u'] % and optional plots can be generated [dplt;'pP']. % % % [1] S. D. Ewert and T. Dau. Characterizing frequency slectivity for envelope fluctuations. % J. Acoust Soc Amer, 108 (3): 1181?196, Sept. 2000. % [2] M. L. Jepsen, S. D. Ewert, and T. Dau. A computational model of human auditory signal processing and perception. % J. Acoust Soc Amer, 124 (1): 422?38, July 2008. % [3] G. Kim, Y. Lu, Y. Hu, and P. C. Loizou. An algorithm that improves speech intelligibility % in noise for normal-hearing listeners. % J. Acoust Soc Amer, 126 (3): 1486?494, Sept. 2009. doi: 10.1121/1.3184603. % [4] J. Tchorz and B. Kollmeier. Estimation of the signal-to-noise ratio with amplitude % modulation spectrograms. Speech Commun., 38: 1?7, 2002. % [5] J. Tchorz and B. Kollmeier. SNR estimation based on amplitude modulation analysis with % applications to noise suppression. % IEEE Trans Speech Audio Processing, 11 (3): 184?92, May 2003. doi: 10.1109/TSA.2003.811542. % bugs: % (3) should scale the derivatives so they get the correct units (e.g. power per second or per mel) % (5) should take care of the scaling so that FFT size does not affect % results % (4) sort out quefrency scale for graphs % (5) could add an option to draw an algorithm flow diagram % Copyright (C) Mike Brookes 1997 % Version: $Id: modspect.m,v 1.11 2009/11/01 21:09:29 dmb Exp $ % % VOICEBOX is a MATLAB toolbox for speech processing. % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % This program is free software; you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by % the Free Software Foundation; either version 2 of the License, or % (at your option) any later version. % % This program is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. % % You can obtain a copy of the GNU General Public License from % http://www.gnu.org/copyleft/gpl.html or by writing to % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% persistent PP if isempty(PP) % Algorithm constants: the first letter indicates the signal domain as follows: % t=time, f=frequency, p=mel freq, m=modulation, q=mel-modulation, d=mod-cepstral PP.logr=120; % maximum dynamic range before taking log (dB) PP.tagc=0; % 1 = do agc PP.tagt=2; % agc time constant PP.tinc=0.25e-3; % time domain frame increment PP.tovl=16; % time domain overlap factor PP.tpad=30e-3; % time domain minimum FFT length (s). must be > 1/min spacing % of mel filters in Hz PP.trnd=1; % 1 to round up time domain FFT to the next power of 2 PP.twin='m'; % time domain window shape: m=hamming, n=hann PP.fpow='p'; % mel filters act on a=amplitude or p=power or c=complex FFT coefs PP.fmin=40; % minimum mel filterbank frequency PP.fmax=10000; % maximum mel filterbank frequency (0 for Nyquist) PP.fwin='t'; % mel filterbank filter shape: t=triangle, m=hamming, n=hann PP.fbin=0.1; % number of mel filters or target spacing in k-mel PP.ppow='a'; % modulation signal is a=amplitude or p=power or l=log-power PP.pinc=16e-3; % modulation frame increment PP.povl=2; % modulation frame overlap factor PP.pwin='m'; % mel domain time-window shape PP.ppad=200e-3; % mel domain minimum FFT length (s). must be > 1/min spacing % of mod-mel filters in Hz PP.prnd=1; % 1 to round up mel domain FFT to the next power of 2 PP.mpow='p'; % transform mel domain a=amplitude or p=power PP.mwin='t'; % mel-mod filterbank filter shape: t=triangle, m=hamming, n=hann PP.mbin=15; % Number of mel-mod filters PP.mmin=50; % minimum modulation frequency PP.mmax=400; % maximum modulation frequency (0 for Nyquist) PP.qpow='a'; % Use mel-modulation domain a=amplitude or p=power or l=log-power PP.qdcq=0; % 1 = Take DCT in modulation direction PP.qdcp=0; % 1 = Take DCT in frequency direction PP.dzcq=0; % 1 = include 0'th coefficient in modulation direction (always included in derivatives) PP.dzcp=0; % 1 = include 0'th coefficient in frequency direction PP.dncp=30; % number of frequency coefficient (excl 0'th) PP.dncq=15; % number of mel-mod coefficient (excl 0'th) PP.ddet=0; % 1 = include delta-time coefficients PP.ddnt=2; % (length of time filter - 1)/2 PP.ddep=0; % 1 = include delta-freq coefficient PP.ddnp=1; % (length of freq filter - 1)/2 PP.unit=0; % 1 = give frequency axes in mels and mod-mels PP.dplt=0; % plotting bitfield: 1=new fig,2=link axes, 4=base coefs, 8=d/dp, 16=d/dt, % 32=waveform, 64=spectrogram, 128=mel-spectrogram, 256=mod-spectrogram, 512 frequency chart PP.cent=0.2; % centile to check for log plotting PP.cchk=0.2; % fraction to require above the centile to allow linear plot end po=PP; % initialize the parameter structure to the default values switch nargin case 0 % list all fields nn=sort(fieldnames(PP)); cnn=char(nn); fprintf('%d Voicebox parameters:\n',length(nn)); for i=1:length(nn); if ischar(PP.(nn{i})) fmt=' %s = %s\n'; else fmt=' %s = %g\n'; end fprintf(fmt,cnn(i,:),PP.(nn{i})); end return case 1 error('no sample frequency specified'); case 2 p=[]; case 3 if isstruct(m) p=m; m=''; else p=[]; end nf=[]; nq=[]; case 4 if isstruct(nf) p=nf; nf=[]; else p=[]; end nq=[]; case 5 if isstruct(nq) p=nq; nq=[]; else p=[]; end end if isstruct(p) % copy specified fields into po structure nn=fieldnames(p); for i=1:length(nn) if isfield(po,nn{i}) po.(nn{i})=p.(nn{i}); end end end % now sort out the mode flags for i=1:length(m) switch m(i) case 'g' po.tagc=1; case {'r','n','m'} po.twin=m(i); po.pwin=m(i); case {'T','N','M'} po.fwin=lower(m(i)); po.mwin=po.fwin; case {'a','c'} po.fpow=m(i); case {'p','l'} po.ppow=m(i); case 'A' po.mpow=lower(m(i)); case {'P','L'} po.qpow=lower(m(i)); case 'f' po.qdcq=1; case 'F' po.qdcp=1; case 'z' po.dzcq=1; case 'Z' po.dzcp=1; case 'd' po.ddet=1; case 'D' po.ddep=1; case 'u' po.unit=1; case 's' po.dplt=bitor(po.dplt,31+512); case 'S' po.dplt=bitor(po.dplt,1023); end end if ~nargout po.dplt=bitor(po.dplt,4); % just plot base coefficients end if ~isempty(nf) po.dncp=nf; end if ~isempty(nq) po.dncq=nq; end lnf=length(nf); if lnf>=1 po.fbin=nf(1); if lnf>=2 po.fmin=nf(2); if lnf>=3 po.fmax=nf(3); if lnf>=4 po.dncp=nf(4); end end end end lnq=length(nq); if lnq>=1 po.mbin=nq(1); if lnq>=2 po.mmin=nq(2); if lnq>=3 po.mmax=nq(3); if lnq>=4 po.dncq=nq(4); end end end end % finally eliminate potential incomaptibilities po.ddep=po.ddep & (po.qdcp==0); % never differentiate along DCT axis po.dzcq=po.dzcq | (po.qdcq==0); % include all coefficients if no DCT po.dzcp=po.dzcp | (po.qdcp==0); % include all coefficients if no DCT logth=10^(-po.logr/10); ns=length(s); axhandle=[]; % % first apply AGC % if po.tagc tau=po.tagt*fs; % time constant for power filtering ax2=[1 -exp(-1/tau)]; bx2=sum(ax2); sm=mean(s(1:min(ns,round(tau))).^2); % initialize filter using average of 1 tau if sm>0 s=s./sqrt(filter(bx2,ax2,s.^2,-ax2(2)*sm)); % normalize to RMS=1 end end if bitand(po.dplt,32) if bitand(po.dplt,1) figure; end plot((1:ns)/fs,s); axhandle(end+1)=gca; title('Input after AGC'); xlabel('Time (s)'); end % % % calculate power spectra % ni=ceil(po.tinc*fs); % frame increment in samples nw=ni*po.tovl; % window length nf=ceil(max(nw,po.tpad*fs)); % FFT length: pad with zeros if po.trnd nf=2^nextpow2(nf); % round up FFT length to next power of 2 else nf=round(nf); end switch po.twin case 'm' w=hamming(nw+1)'; w(end)=[]; % Hamming window ([2] uses Hanning) case 'n' w=hanning(nw-1)'; % Hanning window: underlying period is nw case 'r' w=ones(nw,1); % rctangular window end fx=rfft(enframe(s,w,ni),nf,2); [nft,nff]=size(fx); % number of frames and frequency bins if bitand(po.dplt,64) if bitand(po.dplt,1) figure; end fp=fx.*conj(fx); imagesc(((0:nft-1)*ni+(nw+1)/2)/fs,(0:nff-1)*fs*0.001/nf,10*log10(max(fp,max(fp(:))*1e-6))'); axhandle(end+1)=gca; hanc= colorbar; set(get(hanc,'ylabel'),'String','dB'); axis('xy'); title('Input Spectrogram'); xlabel('Time (s)'); ylabel('Frequency (kHz)'); end [mbk,ffm]=melbankm(po.fbin,nf,fs,po.fmin/fs,min(po.fmax/fs+(po.fmax==0),0.5),po.fwin); switch [po.fpow po.ppow] case 'aa' px=abs(fx)*mbk'; % amplitude domain filtering case {'ap','al'} px=(abs(fx)*mbk').^2; % amplitude domain filtering case 'pa' px=sqrt((fx.*conj(fx))*mbk'); % power domain filtering case {'pp','pl'} px=(fx.*conj(fx))*mbk'; % power domain filtering case 'ca' px=abs(fx*mbk'); % complex domain filtering case {'cp','cl'} px=fx*mbk'; % complex domain filtering px=px.*conj(px); % convert to power end if po.ppow=='l' px=log(max((px),max(px(:))*logth)); % clip before log end if bitand(po.dplt,128) if bitand(po.dplt,1) figure; end switch po.ppow case 'a' imagesc(((0:nft-1)*ni+(nw+1)/2)/fs,ffm/1000,20*log10(max(px,max(px(:))*1e-3))'); case 'p' imagesc(((0:nft-1)*ni+(nw+1)/2)/fs,ffm/1000,10*log10(max(px,max(px(:))*1e-6))'); case 'l' imagesc(((0:nft-1)*ni+(nw+1)/2)/fs,ffm/1000,max(px*10/log(10),max(px(:))*10/log(10)-60)'); end axhandle(end+1)=gca; hanc= colorbar; set(get(hanc,'ylabel'),'String','dB'); axis('xy'); title('Mel Spectrogram'); xlabel('Time (s)'); ylabel('Frequency (k-mel)'); end npf=length(ffm); % number of mel filters mni=ceil(po.pinc*fs/ni); % frame increment in spectral frames mnw=mni*po.povl; % window length mnf=ceil(max(mnw,po.ppad*fs/ni)); % FFT length: pad with zeros if po.prnd mnf=2^nextpow2(mnf); % round up FFT length to next power of 2 else mnf=round(mnf); end nmt=fix((nft-mnw+mni)/mni); % number of modulation spectrum frames ix=repmat((1:mnw)',1,nmt)+repmat((0:nmt-1)*mni,mnw,1); % time index for all modumation spectrum frames mx=rfft(reshape(px(ix(:),:),mnw,nmt*npf),mnf,1); % find modulation spectrum % nmm=size(mx,1); % number of modulation spectrum bins [not needed] [qbk,qqm]=melbankm(po.mbin,mnf,100*fs/ni,po.mmin*ni/fs,min(po.mmax*ni/fs+(po.mmax==0),0.5),po.mwin); % multiply frq by 100 to make it alsmost log nqq=length(qqm); switch po.mpow case 'a' qx=reshape(qbk*abs(mx),nqq,nmt,npf); case 'p' qx=reshape(qbk*(mx.*conj(mx)),nqq,nmt,npf); end switch [po.mpow po.qpow] case 'aa' qx=reshape(qbk*abs(mx),nqq,nmt,npf); % amplitude domain filtering case {'ap','al'} qx=reshape(qbk*abs(mx),nqq,nmt,npf).^2; % amplitude domain filtering + power out case 'pa' qx=sqrt(reshape(qbk*(mx.*conj(mx)),nqq,nmt,npf)); % power domain filtering + amp out case {'pp','pl'} qx=reshape(qbk*(mx.*conj(mx)),nqq,nmt,npf); % power domain filtering end if po.qpow=='l' qx=log(max((qx),max(qx(:))*logth)); % clip before log end tt=((1+mnw*ni+nw-ni)/2+(0:nmt-1)*mni*ni)/fs; % time axis of output frames if bitand(po.dplt,256) && (~bitand(po.dplt,4) || po.qdcq>0 || po.qdcp>0) if bitand(po.dplt,1) figure; end switch po.qpow case 'a' dqx=20*log10(max(qx,max(qx(:))*1e-3)); case 'p' dqx=10*log10(max(qx,max(qx(:))*1e-6)); case 'l' dqx=max(qx*10/log(10),max(qx(:))*10/log(10)-60); end dqx(end+1,:,:)=max(dqx(:)); % insert a border dqx=reshape(permute(dqx,[2,1,3]),nmt,(nqq+1)*npf); ffq=ffm(1)+((0:(nqq+1)*npf)-(nqq-1)/2)*(ffm(2)-ffm(1))/(nqq+1); % mel frequencies imagesc(tt,ffq/1000,dqx'); axhandle(end+1)=gca; hanc= colorbar; set(get(hanc,'ylabel'),'String','dB'); axis('xy'); title('Modulation Spectrogram'); xlabel('Time (s)'); ylabel('Frequency (k-mel)'); end ndq=nqq; % number of coefficients to use in the q direction if po.qdcq dx=reshape(rdct(reshape(qx,nqq,nmt*npf)),nqq,nmt,npf); % take dct in q direction ndq=min(ndq,po.dncq+1); % "+1" to include the 0'th coefficient else dx=qx; end ndf=npf; % number of coefficients to use in the p direction if po.qdcp dx=permute(reshape(rdct(reshape(permute(qx,[3,1,2]),npf,nqq*nmt)),npf,nqq,nmt),[2 3 1]); % take dct in p direction ndf=min(ndf,po.dncp+1); % "+1" to include the 0'th coefficient elseif po.ddep % calculate the frequency derivative nv=po.ddnp; vv=(-nv:nv)*-3/(nv*(nv+1)*(2*nv+1)*(ffm(2)-ffm(1))); dxdp=filter(vv,1,dx,[],3); dxdp=dxdp(:,:,[repmat(2*nv+1,1,nv) 2*nv+1:npf repmat(npf,1,nv)]); % replicate filter outputs at ends end if po.ddet % calculate the time derivative nv=po.ddnt; vv=(-nv:nv)*-3/(nv*(nv+1)*(2*nv+1)*(tt(2)-tt(1))); dxdt=filter(vv,1,dx,[],2); dxdt=dxdt(:,[repmat(2*nv+1,1,nv) 2*nv+1:nmt repmat(nmt,1,nv)],:); % replicate filter outputs at ends end nqj=ndq-(po.dzcq==0); nqk=nqj+ndq*(po.ddep+po.ddet); npk=ndf-(po.dzcp==0); c=zeros(nqk,npk,nmt); c(1:nqj,:,:)=permute(dx(1+ndq-nqj:ndq,:,1+ndf-npk:ndf),[1,3,2]); % base coefficients if bitand(po.dplt,512) if bitand(po.dplt,1) figure; end ffx=repmat(mel2frq(ffm(:)),1,nqq)/1000; qqx=repmat(mel2frq(qqm(:)')/100,npf,1); plot(qqx,ffx,'xb'); axis([[qqx(1) qqx(end)]*[1.1 -0.1; -0.1 1.1] [ffx(1) ffx(end)]*[1.1 -0.1; -0.1 1.1]]); title('Frequency bin centres'); xlabel(sprintf('%.1f : %.1f : %.1f mod-mel = Modulation Frequency (Hz)',qqm(1)/100,(qqm(2)-qqm(1))/100,qqm(end)/100)); ylabel(sprintf('%.0f : %.0f : %.0f mel = Frequency (kHz)',ffm(1),ffm(2)-ffm(1),ffm(end))); end if bitand(po.dplt,4) if bitand(po.dplt,1) figure; end dqx=dx(1+ndq-nqj:ndq,:,1+ndf-npk:ndf); dqxmx=max(abs(dqx(:))); dqxge=dqx>=0; dbfact=2-(po.qpow=='p'); % 2=amplitude, 1=power dqxa=max(abs(dqx),dqxmx*10^(-6/dbfact)); % clip at -60 dB dqxmn=min(dqxa(:)); dblab=''; if(mean(dqxa(:)>dqxmn+po.cent*(dqxmx-dqxmn)))<po.cchk dboff=abs(round(dbfact*10*log10(dqxmn))); dblab='{\pm}dB'; if ~all(dqxge(:)) % not all positive dqxa=dqxa/dqxmn; % force log to be always positive if dqxmn>1 && dboff~=0 dblab=sprintf('{\\pm}dB - %d',dboff); else dblab=sprintf('{\\pm}dB + %d',dboff); end else dblab='dB'; end dqx=dbfact*10*log10(dqxa).*(2*dqxge-1); end dqx(end+1,:,:)=max(dqx(:)); % insert a border dqx=reshape(permute(dqx,[2,1,3]),nmt,(nqj+1)*npk); ffq=ffm(1)+((0:(nqj+1)*npf)-(nqj-1)/2)*(ffm(2)-ffm(1))/(nqj+1); % mel frequencies imagesc(tt,ffq/1000,dqx'); axhandle(end+1)=gca; hanc= colorbar; set(get(hanc,'ylabel'),'String',dblab); axis('xy'); if po.qdcq title('Modulation spectrum DCT'); else title('Modulation spectrum'); end xlabel('Time (s)'); if po.qdcp ylabel('Quefrency (k-mel^{-1})'); else ylabel('Frequency (k-mel)'); end end if po.ddep c(nqj+1:nqj+ndq,:,:)=permute(dxdp(1:ndq,:,1+ndf-npk:ndf),[1,3,2]); % add on p delta if bitand(po.dplt,8) if bitand(po.dplt,1) figure; end dqx=dxdp(1:ndq,:,1+ndf-npk:ndf); dqxmx=max(abs(dqx(:))); dqxge=dqx>=0; dbfact=2-(po.qpow=='p'); % 2=amplitude, 1=power dqxa=max(abs(dqx),dqxmx*10^(-6/dbfact)); % clip at -60 dB dqxmn=min(dqxa(:)); dblab=''; if(mean(dqxa(:)>dqxmn+po.cent*(dqxmx-dqxmn)))<po.cchk dboff=abs(round(dbfact*10*log10(dqxmn))); dblab='{\pm}dB'; if ~all(dqxge(:)) % not all positive dqxa=dqxa/dqxmn; % force log to be always positive if dqxmn>1 && dboff~=0 dblab=sprintf('{\\pm}dB - %d',dboff); else dblab=sprintf('{\\pm}dB + %d',dboff); end else dblab='dB'; end dqx=dbfact*10*log10(dqxa).*(2*dqxge-1); end dqx(end+1,:,:)=max(dqx(:)); % insert a border dqx=reshape(permute(dqx,[2,1,3]),nmt,(ndq+1)*npk); ffq=ffm(1)+((0:(ndq+1)*npf)-(ndq-1)/2)*(ffm(2)-ffm(1))/(ndq+1); % mel frequencies imagesc(tt,ffq/1000,dqx'); axhandle(end+1)=gca; hanc= colorbar; set(get(hanc,'ylabel'),'String',dblab); axis('xy'); if po.qdcq title('Modulation spectrum DCT = freq derivative'); else title('Modulation spectrum = freq derivative'); end xlabel('Time (s)'); if po.qdcp ylabel('Quefrency (k-mel^{-1})'); else ylabel('Frequency (k-mel)'); end end end if po.ddet c(nqk-ndq+1:nqk,:,:)=permute(dxdt(1:ndq,:,1+ndf-npk:ndf),[1,3,2]); % add on t delta if bitand(po.dplt,16) if bitand(po.dplt,1) figure; end dqx=dxdt(1:ndq,:,1+ndf-npk:ndf); dqxmx=max(abs(dqx(:))); dqxge=dqx>=0; dbfact=2-(po.qpow=='p'); % 2=amplitude, 1=power dqxa=max(abs(dqx),dqxmx*10^(-6/dbfact)); % clip at -60 dB dqxmn=min(dqxa(:)); dblab=''; if(mean(dqxa(:)>dqxmn+po.cent*(dqxmx-dqxmn)))<po.cchk dboff=abs(round(dbfact*10*log10(dqxmn))); dblab='{\pm}dB'; if ~all(dqxge(:)) % not all positive dqxa=dqxa/dqxmn; % force log to be always positive if dqxmn>1 && dboff~=0 dblab=sprintf('{\\pm}dB - %d',dboff); else dblab=sprintf('{\\pm}dB + %d',dboff); end else dblab='dB'; end dqx=dbfact*10*log10(dqxa).*(2*dqxge-1); end dqx(end+1,:,:)=max(dqx(:)); % insert a border dqx=reshape(permute(dqx,[2,1,3]),nmt,(ndq+1)*npk); ffq=ffm(1)+((0:(nqj+1)*npf)-(nqj-1)/2)*(ffm(2)-ffm(1))/(nqj+1); % mel frequencies imagesc(tt,ffq/1000,dqx'); axhandle(end+1)=gca; hanc= colorbar; set(get(hanc,'ylabel'),'String',dblab); axis('xy'); if po.qdcq title('Modulation spectrum DCT = time derivative'); else title('Modulation spectrum = time derivative'); end xlabel('Time (s)'); if po.qdcp ylabel('Quefrency (k-mel^{-1})'); else ylabel('Frequency (k-mel)'); end end end if po.unit % frequency units are in mels ff=ffm; qq=qqm/100; else ff=mel2frq(ffm); qq=mel2frq(qqm)/100; end if length(axhandle)>1 if ~bitand(po.dplt,32+64) linkaxes(axhandle) else linkaxes(axhandle,'x') end end